The one-dimensional Bose-Hubbard Model with nearest-neighbor interaction 
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We study the one-dimensional Bose-Hubbard model using the Density-Matrix Renormalization 
Group (DMRG). For the cases of on-site interactions and additional nearest-neighbor interactions 
the phase boundaries of the Mott-insulators and charge density wave phases are determined. We 
find a direct phase transition between the charge density wave phase and the superfluid phase, and 
no supersolid or normal phases. In the presence of nearest-neighbor interaction the charge density 
wave phase is completely surrounded by a region in which the effective interactions in the superfluid 
phase are repulsive. It is known from Luttinger liquid theory that a single impurity causes the 
system to be insulating if the effective interactions are repulsive, and that an even bigger region of 
the superfluid phase is driven into a Bose-glass phase by any finite quenched disorder. We determine 
the boundaries of both regions in the phase diagram. The ac- conductivity in the superfluid phase 
in the attractive and the repulsive region is calculated, and a big superfluid stiffness is found in the 
attractive as well as the repulsive region. 



I. INTRODUCTION 

At zero temperature superfluid-to-insulator |43hase 
transitions can be found in bosonic lattice systemil. Ex- 
perimental realizations of such systems are superconduct- 
ing islands or grains connected by Josephson junctions, 
in which the relevant particles are the bosonic Cooper 
pairs. The transitions at zero temperature belong to the 
class of quantum phase transitions that are not driven 
by thermal, but by quantum fluctuations. These quan- 
tum fluctuations are controlled by system parameters like 
the charging energy of the superconducting islands and 
the Josephson coupling between them. Depending on 
such parameters, the system can assume different forms 
of long-range order. In two dimensions superfluid-to- 
insulator phase transitions airzcro temperature were ob- 
served in thin granular filmgacl and in fabricated Joseph- 
son junction arraysQ'u. 

Recently experiments were carried out in one dimen- 
sional systems of fabricated Josephson junctions. In 
chains one junction wide and 63, 127 and 255 junctions 
long with a tunable Josephson coupling a superfluid- 
insulator transition was observedQ. In another experi- 
ment with fabricated Josephson junctions, long and nar- 
row arrays formed an effectively one-diroHnsional lattice 
for lattice fluxes formed by Cooper pairsBcl. The density 
of these fluxes was controlled by an external magnetic 
field perpendicular to the array, and for small ratios of 
Josephson coupling to charging energy, insulating charge 
density wave phases with densities 1/3, 1/2, 2/3, 1.. were 
found. 

In all of these systems the relevant particles, the 
Cooper pairs or the lattice fluxes, are, at least approxi- 
mately, bosonic. In this paper we will study the phase di- 
agram of strongly interacting bosons on one-dimensional 
lattices. In both of the above experiments on one dimen- 
sional Josephson junction arrays the range of the interac- 



tions were reported to be several sites long. To study the 
effect of longer-ranged interactions we will take nearest- 
neighbor interactions into account. 

In the presence of on-site interactions only, Mott- 
insulators are fomid at integer densities, surrounded by a 
superfluid phaseEJ. Additional nearest-neighbor interac- 
tion leads to charge density wave phases at half-integer 
densities. At the transition from the charge density wave 
phase to the superfluid phase two forms of long range 
order are involved: charge density wave order and super- 
fluid order. If there is a direct phase transition from the 
charge density wave phase to the superfluid phase, one 
type of order appears at the same point where the other 
vanishes. In two and higher dimensions an intermediate 
phase, called the supersolid phase, that has both types of 
order, was found between the charge density wpjejahase 
and the superfluid phase in theoretical modelsE2fEil. For 
experimentaLaystcms the existence of supersolids is still 
controversially. In one-dinicnsional|-hosonic models su- 
persolids have not been found so farE3. 

The existence of a normal or metallic phase that has 
neither superfluid nor chacge density wave order was 
claimed in oneEj and twcO dimensional bosonic mod- 
els in the high density limit. While normal phases at 
zero temperature have not been ruled out rigorously, it 
has been argued that they should not exisillj. We will 
determine whether normal or supersolid phases exist in 
the one dimensional case. 

Effects of disorder and impurities can be especially 
strong in one-dimension. If the effective interactions are 
repulsive, a single impurity_can make some regions of the 
superfluid phase insulatinglZI, while an even larger area of 
the superfluid phase is titmcd into a lattice glass by any 
finite quenched disordeill3EI. The regions in wiiichjpipu- 
rities and disorder become relevant are knownliS'trEZI, and 
can be determined in the phase diagram without actually 
adding impurities or disorder to the model. Since disor- 
der and impurities are important in experiments, we will 



determine these regions in the phase diagram. We will 
also compare the conductivity and Drude weight for the 
attractive and repulsive regions of the superfluid phase 
in pure systems without impurities or disorder. 

The outline of this paper is as follows: in Section O the 
basic phase diagram and the possible phase transitions of 
the Bose-Hubbard model are discussed. Some aspects of 
the density-matrix renormalization group are discussed 
in Section. Ill . The calculation of the phase boundaries 
is presented in Section IV. The correlation functions in 



the different phases are shown in Section M. In Section VI 
the phase diagram with on-site interaction is presented. 
The possible existenc e of normal or supersolid phases is 
discussed in Section VII, and the phase diagram with 
nearest- neighbor interaction is determined. In Section 



VIII the ac-conductivity and the superfluid stiffness is 



calculated for the Mott-insulator and different regions of 
the superfluid phase. Conclusions are given in Section 

ixl. 



II. THE BOSE-HUBBARD MODEL 

The basic physics of interacting bosons on a lattice 
is contained in the Bose-Hubbard modeU. We use an 
extended version which includes nearest-neighbor repul- 
sion: 
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Here Eq is the energy of the insulator groundstate, E^ 
is the energy of a state with the density of the ground- 
state and an additional particle and E^ of that with an 
additional hole. These energies can be calculated using 
DMRG, which will be discussed in Section IV. Note that 



the chemical potentials /z^ and /x^ are not equal to each 
other, the compressibility of the insulator is k = dp/dfj, 
is zero. The superfluid phase is compressible, and for 
states in the superfluid phase uS = n^. The values of ^p 
and ^^ at t = shown in Fig. |l| can be easily calculated 
analytically. 

At the phase transitions from the insulator to the su- 
perfluid phase where the density remains an integer, the 
model is in the universality iCla^ of the xy-model, and 
there is a Kosterlitz-ThoulessEjtJ phase transition. This 
transition is purely driven by phase fluctuations that are 
determined by t. The particle-hole excitation gap at the 
Kosterlitz-Thouless transition closes as: 
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giving the insulating regions a very pointed shape. The 
commensurate-incommensurate phase boundaries can be 
determined directly by calculating the particle and hole 
excitation energies (Eq. (g) and Eq. (||)), and in principle 
the Kosterlitz-Thouless transition could also be found by 
locating the t at which Eg is zero. But since the energy 
gap closes very slowly (Eq. (0)), small errors in the en- 
ergies lead to a big error in the location tc of the critical 
point. Instead, we will study the decay of the correlation 
functions to find the critical point. 



where the bi are the annihilation operators of bosons on 
site i, rii — bjb^ is the number of particles on site i, and t 
is the hopping matrix element. U is the on-site Coulomb 
repulsion and V is the nearest-neighbor repulsion. The 
energy scale is set by choosing U = I. 

For small interactions or large t the bosons are com- 
pletely delocalized, the system is in a superfluid phase. If 
the density is commensurate with the lattice, and there 
is an interaction with the corresponding wavevector, the 
bosons become localized at small t. In the presence of 
on-site interaction only (T^ = 0), Mott-insulating regions 
with integer density are found. Fig. || is a sketch of 
the phase diagram showing Mott-insulating regions sur- 
rounded by the superfluid phase. 

On most of the phase boundaries between the insu- 
lating phases and the superfluid phase the density of the 
system changes as the phase boundary is crossed from the 
incompressible insulator to the compressible superfluid. 
The location of this commensurate to incommensurate 
density transition can be directly determined as the en- 
ergy it takes to add a particle or hole to the insulator: 
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FIG. 1. Schematic phase diagram with V — 0. It shows 
the Mott-insulators (MI) with density p — 1 and p — 2, sur- 
rounded by the superfluid phase(SF). The Luttinger Liquid 
parameter K is shown at the commensurate-incommensurate 
phase transitions and the Kosterlitz-Thouless transitions at 
the tips of the insulating regions. 



The superfluid phase of interacting bosons in one 
dimension has a hnear dispersion relation for smaU 
wavevectors q and no excitation gap. The low eaeEgM 
physics of this phase is that of a Luttinger LiquidElo'Ell 
with the basic Hamiltonian 



Ho = -— I dx 
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Here n(a;) are density fluctuations and <i>(x) phase fluc- 
tuations (b'^ix) — \/p{x) e'*^"^^), V is the second sound 
velocity rand K determines the decay of the correlation 
functioned: 
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for r ^ p. The interactions and the lattice introduce an 
extra term: 



Hr, 



dxcos{2n<d{x) + 2Trnx{p — po)) , (8) 



where dxQ{x) = tt[po + n(a;)], p is the density of the 
system, po is the density of the insulator (e.g. pa = 1 
for the Mott insulator), and n is the denominator of 
the density of the insulator: po = m/n. For K > Kc 
this term becomes relevant and drives the system into 
an insulating phase. At the Kosterlitz-Thouless tran- 
sition Kc = \ and K^ ^ Mr at the commensurate- 



incommensurate transitiontJ 

At the Mott-insulators with integer densities, the de- 
nominator of the density is n = 1. At the sides of 
the insulator the Luttinger-Liquid parameter is isT = 1, 
and at the Kosterlitz-Thouless transition at the tip it is 
K = 1/2. The parameter K can be determined from 
the correlation functions, and we will locate the tc of the 
Kosterlitz-Thouless transition by finding the t at which 
K = 1/2. 

If additional nearest-neighbor interactions are included 
in the model, charge density wave phases are found 
at half integer densities. They also have a Kosterlitz- 
Thouless transition at the tips, giving them a similar 
shape as the Mott-insulators (Fig. 0). Since the denom- 
inator of their density is n = 2, the parameter K — A 
at the sides of the phase boundary of the charge density 
wave, and A' = 2 at the Kosterlitz-Thouless transition at 
the tip. 

The possible existence of an intermediate phase, super- 
solid or normal, between the charge density wave phase 
and the superfluid phase will be addressed in Section VI] . 



III. DMRG 

To determine the energies and correlation func- 
tions we p^e-the Density-Matrix Renormalization Group 
l3, a numerical method capable of delivering 
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precise results for groundstate properties of low dimen- 
sional strongly interacting system. We use the finite-size 
version of the DMRG algorithm, in which the system is 
built up to a certain size, and the basis of the system is 
then optimized to represent the chosen target states by 
sweeping through the system repeatedly until the basis 
is converged. 

The density matrix weight of the states discarded in a 
DMRG step is a measure of the numerical errors caused 
by the truncation. We found this truncation error to 
depend on the correlation length in the system. At a 
fixed number of states kept, we find very small trunca- 
tion errors in the insulating phases, that grow as the 
phase transition to the superfiuid phase is approached, 
and are biggest in the superfluid phase. Note that in one 
dimension the whole superfiuid phase is critical with a 
diverging correlation length, but the correlation length is 
always finite in finite systems. 

In each DMRG calculation for a given set of model 
parameters we first use the groundstate, the state with an 
additional particle and the state with an additional hole 
as target states. To obtain adequate numerical accuracy 
in all cases, we require the density matrix weight of the 
truncated states A < 5 x 10~^ (see Appendix H). The 
energies of these states are used to calculate the chemical 
potentials (Eq. (|) and Eq. (|)). 

For further sweeps only the groundstate is used as a 
target state. We require the weight of the truncated 
states A < 10~^, and the number of states kept is in- 
creased if necessary. After the basis is converged, which 
usually takes two sweeps, the groundstate correlation 
functions are calculated. 

At the same parameters and number of states, the 
truncation error in a system with periodic boundary con- 
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FIG. 2. Schematic phase diagram with V = 0.4. It shows 
the Mott-insulator(MI) with density p = 1, the charge density 
wave (CDW) phases at densities p — 1/2 and p = 3/2, and 
the surrounding superfluid phase(SF). The K are Luttinger 
Liquid parameter at the phase transitions. 



ditions is usually much higher than with open boundary 
conditions, therefore we use open boundary conditions. 
To keep boundary effects small we add additional terms 
on the boundaries to the Hamiltonian: 



Hboundary ^ V p fli + V p flL] ■ 



(9) 



With this additional term, a particle on the boundary on 
average has the same potential energy as in the rest of 
the system. 



IV. PHASE BOUNDARIES 

As pointed out in Section O, the phase boundaries with 
the exception of the Kosterlitz-Thouless transition can be 
determined by calculating the particle and hole excitation 
energies (Eq. (g) and (H)). Using DMRG we calculate 
these energies in finite systems. We expect quadratic 
system size dependence of the energies of the insulator 
groundstates, and linear system size dependence for the 
states with additional particles and holes. Fig. ^(a) and 
Fig. Q show that the leading term in the scaling of /i^ and 
u^ in the insulator and the superfluid phases is 1/L. Fig. 
y(b) shows pP of the Mott-insulator at p = 1 with V = 
and i = 0.1. In this case the system size dependence is 
very weak, and on the scale of Fig. ||(b) the quadratic 
part of the scaling can be seen. Since the quadratic and 
higher parts contribute only very weakly to the scaling, 
we ignore them and use linear extrapolations from the 
finite system sizes to determine p^ and p^ in the ther- 
modynamic limit. In the insulator phase p^ y^ p^^ since 
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FIG. 3. a) System size dependence of the chemical poten- 
tial in the Mott insulator ("*") at p = 1 with i = 0.1 and 
V = and the charge density wave phase ("o") at p = 1/2 
with i = 0.1 and V = 0.4. The upper set of data points (de- 
noted by p) corresponds to the energy necessary to add an 
additional particle (pc), the lower one (denoted by h) to that 
of adding a hole (/ij?)- The dashed lines show linear fits to 
the data, b) p'' for p = 1, V = and i = 0.1 on an expanded 
scale. 
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FIG. 4. System size dependence of the chemical potential 
in the superfluid phase at f = 0.1. As in Fig. H, the upper 
sot of data points (p) corresponds to the energy necessary to 
add an additional particle (/x?), the lower one (h) to that of 
adding a hole (/ic). The dashed lines are linear fits. 



there is a finite gap Eg — pP. — p^ (Eq. (^)). In the 
superfluid phase the extrapolations for p^ and p'^ should 
result in the same p since Eg = and the system is com- 
pressible. In Fig. H small deviations from this can be 
seen, which we will ignore in our analysis. 



V. THE CORRELATION FUNCTIONS 

A. The local density 

Since open boundary conditions are used, special care 
has to be taken to reduce boundary effects. The most ob- 
vious form of these are local density oscillations. In the 
superfluid phase they show the power-law decay away 
from the edge of the system characteristic for the Lut- 
tinger Liquid. If the density of the bosons is given as a 
rational number p = n/m, the wavelength of the oscil- 
lations is the denominator m of the density - the same 
wavelength as in the density-density correlation functions 
(Eq. (R)). Fig. S shows the local density in a system 
at density p — 3/4 with only on-site interactions, and 
p — 1/4 with additional nearest-neighbor interactions, 
both in the superfluid phase. 

In the insulators the boundary effects decay expo- 
nentially. In addition to the boundary induced den- 
sity oscillations, in the charge density wave phase at 
density p — 1/2 (Fig. ||) there is real long range or- 
der in (rii) in the form of a charge density wave with 
(rii) = p + 5^(— 1)', where 5^ is the structure factor. In 
an infinite system the long range order is due to spon- 
taneous symmetry breaking - in finite systems with even 
numbers of sites we allow this by adding a small symme- 
try breaking term to the chemical potential on the left 
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FIG. 5. Local density (n,;) in the superfluid phase. The 
systems are L = 256 sites long. 



boundary. Without this symmetry breaking term reflec- 
tion symmetry would cause the groundstatc to be the lin- 
ear combination of two charge density wave phases with 
a phase difference of tt, canceling out the long ranged 
oscillations in (n^). Breaking the symmetry in the finite 
system reduces the Hilbert space necessary to represent 
the groundstate by half and also leads to a better conver- 
gence of the DMRG. In the superfluid phase the symme- 
try breaking term just modifies the density oscillations 
at the boundary. 
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FIG. 7. {bi-bo) correlation function in superfluid systems 
with density p = 1/4 at t = 0.1, U = 1 and V = 0.4 and differ- 
ent system sizes L. As the system size is increased the bound- 
ary effects become weaker, and the curves for different system 
sizes look similar over increasing regions. The inset shows the 
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decay parameter K from power-law fits (blbo) ~ r 
16 < r < 32 for the different system sizes (L = 32, 64, 96, 128 
and 256). As the system size grows, K converges to the value 
for infinite systems. The dashed lines show the extrapolation 
to the upper and lower limit of K. The mean value of the up- 
per and lower limit is taken as K, and the difference gives an 
estimate of the error. In the shown case it is A" = 1.40 ±0.03. 



B. The hopping correlation function 
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FIG. 6. Local density fiuctuations in the charge density 
wave phase with density p = 0.5 at t — 0.1, U = 1 and 
V = 0.4 in a L = 256 system. The fluctuations induced by 
the right boundary at i = 256 decay quickly, but the true long 
range oscillations go throughout the system. 



In the superfluid phase the hopping correlation func- 
tion r(r) = Cb^j.bo) decays with the power-law behavior 
given in Eq. (^, which can be used to determine the Lut- 
tinger Liquid parameter K. As discussed above, there are 
local density fluctuations (n^) in the finite systems. The 
creation operator can be represented by a density (n.i) 
and a phase (j)i part: b] = \/ {rii) exp {i(j)i). As discussed 
above, there are local density fluctuations {rii) in the fi- 
nite systems, and they will affect the correlation function 
{b\b-). Since the local density oscillations, with the ex- 
ception of the charge density wave phase, are boundary 
induced, and we are interested in the properties of an 
infinite system, we reduce the effect of the local density 
oscillations averaging over pairs of {b\b,) with | i — j |= r. 
To minimize boundary effects we place i and j symmet- 
rically around the center. 

Fig. |7| shows the power-law behavior in the superfluid 
phase for small r, which is modified to a faster decay 
closer to the system boundaries. The bigger the systems 
are, the bigger is the region in which the correlation func- 
tions fit the algebraic decay, and also the region in which 
the correlation functions look the same for different sys- 
tem sizes. In all cases the two biggest systems we calcu- 
late are at least 128 and 256 sites long. To estimate K , 
we fit a-r~^l'^ to the numerical data for 16 < r < 32. For 



systems with 128 and 256 sites the boundary effects in 
this region are small, while the distance is also big enough 
to avoid short-ranged (non-Luttinger Liquid) effects. 

With increasing system size the boundary effects get 
weaker, resulting in a decreasing K that asymptotically 
approaches the infinite size value. To find a simple es- 
timate of this we use the K determined in the biggest 
system as an upper limit K^, and the linear extrapola- 
tion from the values in the two biggest systems as a lower 
limit Ki. We take the mean value K = (Ku + Ki)/2 and 
estimate the error as AK = {K^ ~ Ki)l2. 



C. The density-density correlation function 



VI. ON-SITE INTERACTIONS 

Using the methods described above we determine the 
phase diagram in the presence of on-site interactions only. 
Fig. shows the Mott-insulator with density p = 1, sur- 
rounded by the superfluid phase. In Fig. [l^ the tip of 
the insulator is shown on an expanded scale. The very 
pointed tip of the insulating region reflects the closing of 
the energy gap given by Eq. (p. Fig. g also shows re- 
sults for the commensurate-incommensurate phase tran- 
sition from twelfth order perturbation theory. The excel- 
lent agreement with DMRG confirms the high accuracy 



The density-density correlation functions are calcu- 
lated in the same way as the hopping correlation function. 
However, in this case it is necessary to subtract the static 
expectation values, measuring (riinj) — (ni){nj), instead 
of just taking {niUj). 

Fig. H shows the density-density correlation function 
at density p — 1/4. A fit with Eq. (M) works fairly well, 
but the first term -^(27rpr)~^ could not be observed. In- 
stead of the correlation function being bigger for small r, 
we find it to be smaller. Eq. (0) only necessarily holds 
at large distances, and the short range behavior we see is 
dominated by the repulsive interaction between the par- 
ticles. In fitting A{pr)^^^-^ cos 2TTpr to the data, a cut-off 
at small distances has to be made. While this works well 
enough to confirm that the correlation functions decay 
with a power-law behavior, the uncertainties in the fit 
are too high to determine K. 
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FIG. 10. Same as |9[ with the tip of the Mott-insulator on 
an expanded scale. 
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FIG. 11. if vs t at the Kosterlitz-Thouless transition 
from the density one Mott-insulator to the superfluid phase. 
The critical point is determined by finding the t at which 
K — 1/2. The biggest system sizes are L — 256 sites except 
for t = 0.29, 0.3, 0.31, where they are L = 1024. 



achieved. 

To find the critical point of the Kosterlitz-Thouless 
transition at the tip of the Mott-insulator, we d eterm ine 
the t at which K = 1/2. As described in Section |VB| , the 
K are determined by fitting power-law behavior to the 
decay of the hopping correlation function r(r) = (6j6o). 
Due to finite-size effects, the result can depend on the 
interval of r that is used for this fit. While we found 
16 < r < 32 to be reasonable in general, there are loga- 
rithmic finite-size effects at the Kosterlitz-Thouless tran- 
sition that require a more detailed inspection. Fig. |ll| 
shows K^s that were determined with different fitting in- 
tervals plotted versus t. 

For the t = 0.29, t ^ 0.3 and t = 0.31 systems with 
L = 512 and L — 1024 sites were calculated to keep 
finite-size effects small. Table H shows the tc found with 
different fitting intervals. Due to the finite-size effects the 
tc go to higher t as the fitting interval is shifted to bigger 
r. While using a fitting interval of 16 < r < 32 is a good 
compromise between avoiding finite-size effects by using 
small r, and finding an asymptotic value for the decay 
by using big r, the error determined by this fit alone 
underestimates the true error. With an error estimated 
from the effect of the choice of the fitting interval, we 

TABLE 1. The location of the critical point tc depending 
on the interval of r used for the fits to r(r) = {blbo). 

r tc 



4 < r < 8 
8 < r < 16 
16 < r < 32 
32 < r < 48 
48 < r < 64 



0.2874 ± 0.0001 
0.2938 ± 0.0001 
0.2968 ± 0.0003 
0.3062 ± 0.0003 
0.3107 ±0.01 



findic = 0.297 ±0.01. 

Determination of the Kosterlitz-Thouless transition 
was attempted in several previous studies. For a trun- 
cated model with a maximum of two particles per site 
the critical point was found at if ^ = l/(2\/3) « 0.289 
with the Bethe-AnsatfJ. In a combination of exact- 
diagonalization for system with up to L = 12 sites and 
renormalization gijoup Kashurnikov and Svistunov found 
tc — 0.304 ± 0.002!^, and together with Kravasin found 
tc = 0.300 ± O.OOSEI in a quantum Monte Carlo study. 
An exact diagonalization approacij-|rcported the critical 
point to be at tc = 0.275 ± O.OOsN, anjd_a study using 
12th order strong coupling expansionfc3 located it at 
tc = 0.26 ±0.01. 

There were two previous studies using DMRG to locate 
the Kosterlitz-Thouless tr|aiisition in the Bose-Hubbard 
model. In the first studyEZI a particle cut-off of n = 4 
per site was used, and the critical point was determined 
by using a phase twist = tt and periodic boundaries to 
calculate the superfluid stiffness. The criticaLpoint was 
found at tc = 0.298. In another DMRG studyB, the crit- 
ical point was found at tc — 0.277 ± 0.01 by determining 
the t at which K — 1/2, similar to the procedure used in 
the present work. In both of these DMRG studies the nu- 
merical accuracy was limited due to the use of the infinite 
size version of the DMRG algorithm, periodic boundary 
conditions and a small number of states in the DMRG 
basis. The biggest system sizes were smaller than L — 80. 

The range of these results demonstrates the diffi- 
culty involved in determining the critical point of the 
Kosterlitz-Thouless transition, which is mostly due to 
logarithmic finite-size effects close to the critical point. 
The large system sizes used in this paper should com- 
pensate for this within the given error bars, and yield a 
reliable result. 

The phase diagram shown in Fig. has a very in- 
teresting feature. Imagine moving on a line of constant 
chemical potential /i, for example /i = 0.15, and starting 
at small t, moving toward bigger t. The particle density 
along this line is illustrated in Fig. |l^. For small t the 
system is in the Mott-insulator phase. At i « 0.1 there is 
a phase transition to the superfluid phase, with densities 
p < \. The density decreases up to a minimum, then it 
start increasing again. At t « 0.26 the density goes up 
to p = 1 again, and there is another phase transition, 
this time reentering the Mott-insulating phase from the 
superfluid phase. Increasing t further leads to another 
phase transition from the Mott-insulator to the super- 
fluid phase, this time with p > 1, and the density in- 
creasing further with increasing t. 

To gain more insight intcLlhis, we compare to results 
from a mean field approacl£j. Fig. 13 shows the phase 
diagram with the Mott-insulator at density p = I, sur- 
rounded by the superfluid phase. In the superfluid phase, 
the lines of constant density slope downward as t is in- 
creased. This is not only found in one dimension, but in 
all dimensions. The limit oi t —> oo corresponds to keep- 
ing t constant and setting the interactions to zero. If the 
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FIG. 12. Illustration of the phase transitions between the 
Mott-insulator (MI) and the superfiuid phase (SF) on a line 
of constant chemical potential [s, — 0.15. 
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FIG. 14. The density on a line of constant chemical po- 
tential /i = 0.3 in the Mean-Field phase diagram (Fig. |l3|). 



interactions arc zero the system goes from a superfiuid 
phase to a Bose-Einstein condensate, in which every par- 
ticle has an energy of —2t. If the chemical potential is 
smaller than — 2i, the system is empty, because it costs 
energy to put a particle in, and for chemical potentials 
bigger than — 2i the number of particles goes to infinity, 
because every additional particle reduces the total energy 
of the system. Going back to the picture of constant in- 
teractions and changing t, this means that the density of 
the system always goes to infinity as t is increased. 

In dimensions two and higher the superfluid-insulator 
transition on the line of constant density is a second order 
transition, and the tip of the insulating region is round. 
Fig. 14 shows the density on a line of constant chemical 
potential /i = 0.3. At the phase transition from the Mott- 
insulator to the superfiuid phase the density first drops 
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FIG. 13. The Mean-Field phase diagramE3 in dimension- 
less units, showing lines of constant density. 



as t is increased, and then increases again. 

In one dimension the tip of the insulating region is very 
long and narrow due to the Kosterlitz-Thouless transi- 
tion, and it is possible to reenter into the insulator at 
p=l. 



VII. NEAREST-NEIGHBOR INTERACTION 

Longer ranger interactiops have been found to be im- 
portant in experimentsQia. We now include nearest- 
neighbor interactions by setting V — 0.4. Due to the 
nearest-neighbor interactions a new insulator phase ap- 
pears at half integer densities. It is a charge density wave 
phase (CDW) with a wavelength of two sites, and like the 
Mott-insulator at integer density it has an excitation gap 
and is incompressible. The crystalline order is character- 
ized by a non-zero structure factor 



s^-^Ei-^y-'^i^ 



(10) 



In Fig. ^ the local density oscillations in the charge den- 
sity wave phase are shown. A small boundary effect can 
be seen, but the main feature arc long-range density os- 
cillations throughout the system that do not decay. An 
order parameter {rn — p) can be defined to describe this 
charge density wave, even in one dimension. 

In one dimension the superfiuid phase is signalled by 
a diverging correlation length 



r r 

and a non-zero superfiuid stiffness 

Ps = hm L -— , 

<t>^Q dm-' 



(11) 



(12) 



which is proportional to the Drude weight D ^ ps- In two 
and higher dimensions there is also an order parameter 
(6j) 7^ 0. In one dimension the whole superfluid phase is 
critical, and there is no order parameter. 

At the transition from the charge density wave phase 
to the superfluid phase both types of order are involved: 
the crystalline order in the charge density wave phase and 
the superfluid order in the superfluid phase. In addition 
to a direct phase transition from the charge density wave 
to the superfluid at which the crystalline order vanishes 
at the same point where the superfluid order appears, 
there is the possibility of an intermediate phase. Tab. 
H shows the possible phases close to density p = 1/2 in 
a bosonic system with on-site and nearest-neighbor in- 
teraction in two or higher dimensions. In addition to 
the charge density wave and the superfluid phase, su- 
persolids that have both-ferms of order were found p. 
two dimensional modelsEjO. Baltin and Wagenblasttil 
found a region that has neither superfluid stiffness nor 
charge density wave in a one dimensional bosonic model 
in the high density. The possible existence of such a phase 
was also recently predicted for a two dimensional bosonic 
model in the high density limit by Das and Doniachllj, 
who call it a Bose-metal. 

In Appendix O strong coupling expansions are used 
to illustrate the difference between the commensurate- 
incommensurate phase transition at p = 1/2 in one and 
two dimensions. Strong coupling expansions can be used 
to study the insulator, but not the superfluid. To study 
the low-energy behavior of the superfluid phase the Lut- 
tinger liquid can be used. In addition to the basic Lut- 
tinger liquid Hamiltonian (Eq. (^), the lattice and the 
interactions introduce scattering terms (Eq. d|)). These 
only contribute at p = 1/2, where they can drive the sys- 
tem into a different phase, but not at nearby densities. 
At incommensurate densities close to p ~ 1/2, the wave- 
function is incommensurate with the lattice and hence 
cannot be pinned to the lattice to form an insulator. Of 
course this would be changed if there were impurities or 
disorder, but a pure system in one dimension is in the 
Luttinger liquid phase unless it is at a density commen- 
surate with the lattice and the interactions. 

At density p = 1/2 DMRG can be used to determine 
if there is an intermediate phase or a direct phase tran- 
sition from the charge density wave phase to the super- 
fluid phase. To do this, we investigate the relationship 
between the superfluid and crystalline order at the phase 
transition. The onset of superfluidity is signaled by a 

TABLE II. Possible phases and their order parameters 
close to density p — 1/2. 
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FIG. 15. Structure factor S-n versus the inverse correlation 
length ^ for diflerent system sizes. The dashed lines show 
power-law fits. 



diverging correlation length ^ (Eq. (pT|) ), charge density 
order is meassured by 5'^ > (Eq. (|lOf)). 

The model is in the universality class of the xy-model 
at the phase transition on the line of constant density 
at p = 1/2. Due to the Kosterlitz-Thouless transition 
expected on this line, it is difficult to determine exactly 
when the structure factor or the inverse correlation length 
go to zero. But it is possible to study the dependence of 
the structure factor on the correlation length. Fig. Il^ 
shows that for small values the structure factor depends 
on the correlation length by a power-law: 



i{s-,)-^-i[s^ = Q)-^^{s^r 



(13) 



To keep the effect of the boundaries small for the calcu- 
lation of both the structure factor 5^ and the correlation 
length ^, only sites that were at least a quarter of the 
system size away from the boundaries were taken into 
account. 

In Fig. nq the extrapolated inverse correlation length 
l/C('S'7r — 0) at zero structure factor is plotted against the 
inverse system size. The linear fit to the three biggest 
system sizes shows that 1/^(5^ — 0) goes to zero for 
infinite systems. From this and the power-law behavior in 
Fig. nS we conclude that there is a power-law dependence 
of the structure factor on the correlation length, and that 
in infinite systems the correlation length diverges at the 
same point at which the structure factor goes to zero. 
This means that there is a direct phase transition from 
the charge density wave phase to the superfluid phase, 
and no supersolid or normal phase in between. 

The phase boundaries of the charge density wave phase 
can be found in the same way as those of the Mott- 
insulator, and we use the methods used for the on-site 
only interaction case to calculate the phase diagram. 
Fig. n^ shows the phase diagram in the region of the 
p = 1/2 charge density wave phase and the p = \ 
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FIG. 16. The inverse correlation length 1/^ at S,r = 
obtained from Fig. hsi versus the inverse system size. The 
dashed line shows a linear fit to the three biggest system sizes. 



Mott-insulator. Like the shape of the Mott-insulator 
with on-site interaction only, the shapes of the insulat- 
ing regions reflect the Kosterlitz-Thouless transitions at 
the tips. The tips are also bending down, allowing re- 
entrance phase transitions from the superfluid to the in- 
sulating phases. We find the Kosterlitz-Thouless transi- 
tions at i*" = 0.404±0.02 for the Mott-insulator, and at 
iCDW _ 0.125 ±0.003 for the charge density wave phase. 
The critical point at the tip of the charge density wave 
phase had been found in a quantum Monte Carlo study 
at i « 0.1. The accuracy of i^^^ is relatively low because 
K only changes very slowly if t is changed close to this 
transition. This iiad already been observed in an ear- 
lier DMRG studya, where the critical points were found 
at if ^^ = 0.118 ± 0.004 and if ^ « 0.325 ± 0.05. That 
study used the infinite size version of the DMRG, and the 
biggest systems were L = 76 sites long. Since there are 
logarithmic finite-size effects at the Kosterlitz-Thouless 
transitions, the values determined in the present work. 



Mott-insulator 




where systems with up to L = 512 sites are used, are 
much more accurate. 

At the phase boundaries of the charge density wave 
phase the Luttinger Liquid parameter is i^ = 4 except 
for the Kosterlitz-Thouless transition at the tip, where it 
is i^ = 2. The charge density wave phase is surrounded 
by a region where K > \. The effective interactions 
in the Luttinger-Liquid are attractive for K < 1 and 
repulsive for ii' > 1. Kane and Fisher showed that in the 
repulsive region, a single impurity turns the system into 
an insulatoiOnJ. An even larger region of the Luttinger 
liquid is driven into a glfWS phase ii K > 2/3 for any 
finite quenched disordeiOH. 

In the phase diagram these regions {K > 1 and 
K > 2/3) are determined by doing calculations for differ- 
ent densities and t. For example, to determine the t for 
which K{t) = 1 a,t a, given density p, we do calculations 
for different t until we find a pair of ii and ^2 that are 
close to each other, and K{ti) > 1 > K{t2). We then 
determine t{K = 1) by linear interpolation. To find the 
boundaries of the repulsive Luttinger Liquid region, we 
calculate t{K = 1) for various densities. 

The lines with K = 1 and K = 2/3 are shown in Fig. 
O. The repulsive Luttinger Liquid {K > 1) region com- 
pletely surrounds the charge density wave phase. Instead 
of going to t = as the Mott-insulator is approached, 
the K = 1 line ends in the side of the Mott-insulating 
region, where K = I for all of the commensurate- 
incommensurate phase boundary. Fig. UTj illustrates how 
the two lines with K = 1 meet at the phase boundary of 
the Mott-insulator. Although we could not obtain more 
detailed results for densities closer to p = 1, we argue 
that the lines with K > 1 bend towards i = as the 
density gets closer to one, while those with 1/2 < K < 1 
bend towards the tip of the Mott-insulator. 

Lines of constant K with 1/2 < K < 1 are discon- 
tinuous at p = 1, where the system is an insulator for 
K > 1/2, with the tipped shape reflecting the Kosterlitz- 
Thouless behavior. Lines with K < 1/2 are round at 
p — I, and do not reflect the Kosterlitz-Thouless behav- 
ior. Analogous to p = 1, we also observe that the lines of 
constant K are round for iiT < 2 and p = 1/2, where the 
charge density wave phase ends in a Kosterlitz-Thouless 
transition at isT = 2. 



,^ K=l/2 



FIG. 17. Illustration of the lines of constant K in the 
phase diagram. The dotted lines indicate K that are slightly 
smaller or bigger than K = 1. 
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FIG. 18. The phase diagram of the Bose-Hubbard model with nearest-neighbor interaction V = 0.4. The Mott-insulator 
(MI) at density p = 1 and the charge density wave phase (CDW) at density p = 1/2 are surrounded by the superfluid phase. 
The long-dashed lines show the lines of constant density. The solid lines are lines of constant K. The K = 2 line crosses the 
density p = 1/2 line at the Kosterlitz-Thouless (KT) transition at the tip of the charge density wave phaKe. In the region left 
of the K = 1 line, where Tf > 1, the superfluid phase is turned into an Jaajulator by a single impuritytj. In the presence of 
disorder the region left oi K = 2/3 line is turned into a Bose-glass phasetJiJ. The Kosterlitz-Thouless (KT) transition at the 
tip of the Mott-insulator is at the point where the K — 1/2 line line intersects the p = 1 line. 
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VIII. AC-CONDUCTIVITY 

The repulsive region of the Luttinger liquid is turned 
into an insulator by a single impuritjU. This raises the 
question if there is a qualitative difference between the 
conductivity in the repulsive region and the attractive 
region of the Luttinger liquid in the pure system. The 
regular part of the conductivity is given by: 



^-s(^) = 1 ^ 



\{m\j,=o\0)t 



L '^ E„i — Eq 

1 



-Im lira 

luttL 77^0+ 



(0 1 4=0: 



1 



"i='> uj + Eq - H + ir] 
and the current operator is 



■5{U0 - {Era - Eo)) 



.h=o I 0) , (14) 



jq = it 



E' 



-tq7i I 



(&n + l^n - ^-C-) 



(15) 



Recent developments with DMRG make the calcu- 
lation of dynamical |-cprrelation functions like the ac- 
conductivity possibleEa. The conductivity at a frequency 
z = w + iry can be calculated as the direct product of the 
current operator applied to the groundstate 



I jq=o) = jq=0 I 0) , 



and a correction vector 



x{z)) 



1 



Eo- H + iT] 



h=o) 



(16) 



(17) 



By using these two states and the groundstate |0) as 
DMRG target states, the conductivity (j^'^^{z) can be cal- 
culated very precisely. To calculate the conductivity over 
an interval of width rj ranging from wi to uj2, we use cor- 
rection vectors | x{lji -f iri)) and | a;(ti;2 + iv)) ^s target 
states. At the end of the DMRG calculation, when the 
DMRG basis is optimized to represent these states, we 
calculate the conductivity from wi to uj2 ■ Repeating this 
procedure for neighboring intervals, we piece together the 
conductivity for a whole range of frequencies. 

The finite broadening rj in the correction vectors is 
only used to obtain appropriate DMRG target states. To 
calculate the spectrum within the DMRG basis, we use 
a Lanczos method that yields approximate eigenstates 
of the Hamiltonian. The broadening used in our plots is 
then applied to the discrete peaks found with the Lanczos 
method, and is only used for better visualization. 

DMRG calculations work best with open boundary 
conditions. How the current operator is applied in a 
system with open boundary conditions is discussed in 
Appendix O. 

The conductivity in the Mott-insulator phase is shown 
in Fig. 11^. There is an energy gap of Aw = 0.54, with a 
big peak after that and only small excitations at higher 



TABLE in. The Drude weight D for different system 
sizes in the Mott-insulator with density p — 1, t — 0.05, 
U ^ 1,V ^ 0.4. (T) is the kinetic energy, J al''^{uj)duj the 
integral over the ac- conductivity. Also shown is the number 
of states m, the broadening rj of the correction vectors, and 
the truncation error A. 



L 


D/t 


{T)/itL) 


2/f/ar(c. 


)dLU 


V 


m 


A 


32 


0.0084 


0.6392 


0.6307 




0.05 


128 


10-' 


64 


-0.0040 


0.6392 


0.6432 




0.05 


128 


10-" 


128 


-0.0075 


0.6392 


0.6467 




0.05 


128 


10-' 


256 


-0.0036 


0.6392 


0.6428 




0.05 


128 


10-" 



energies. Since it is an insulator phase, we expect to find 
no Drude weight. With the kinetic energy defined as 



{T)=tJ2{bi+ibn + h.c.) 



the Drude weight is given by 
1 



D 



L 



(T) 



duj aY^^ (lu) 



(18) 



(19) 



Note that the Drude weight is proportional to the super- 



m 



fluid stiffness ps given in Eq. 

Since the kinetic energy (T) can be calculated directly 
with DMRG, and we expect D — 0, this is an opportunity 
to verify the consistency of the calculation. In Tab. 
the Drude weight is shown for various system sizes. 



nil 
A 



small finite-size effect can be seen in the data. From the 
differences in the individual values we estimate the error 
of the Drude weight AD = 0.02t, or 2% of ~{T)/L. 

In the superfluid phase we find precursor peaks at small 
frequencies in the conductivity. They are due to the finite 
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FIG. 19. The conductivity a\'^^{ijj) in the Mott-insulator 
at density p = 1, f = 0.05, U = 1, V = 0.4, system size 
L — 256, m — 128 states, two correction vectors as target 
states with rj = 0.05, and broadening rjg — 0.01 for the plot. 
The inset shows the same data on an expanded scale. 
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TABLE IV. The Drude weight D for different system sizes 
in the repulsive region of the superfluid phase at density 
p = 3/4, t = 0.05, U = 1,V = 0.4. Same notation as in 
Tab.fn 



L D/t {T)/{tL) 2/t J al^^Lo)duj r) m 



A 



16 0.776 


0.874 


0.099 


0.2 


128 5 X 10^' 0.2 




32 0.791 


0.885 


0.094 


0.2 


128 7 X 10"^ 0.1 


3 


64 0.799 


0.890 


0.091 


0.2 


128 5 X 10-*^ 0.08 


£ 


64 0.798 


0.890 


0.092 


0.05 


256 1 X 10"^ 0.08 


D 

■7 


128 0.799 


0.892 


0.093 


0.2 


128 7 X 10"^ 0.02 


o 


128 0.801 


0.892 


0.091 


0.1 


256 6 X 10^^ 0.02 




128 0.793 


0.892 


0.099 


0.05 


128 1 X 10"^ 0.02 




128 0.797 


0.892 


0.095 


0.05 


256 1 X lO"** 0.02 




256 0.800 


0.893 


0.093 


0.2 


512 5 X 10'^ 0.02 





widthEj of the wavevector q, which is Aq = 4^/5/ L. Fig. 
EQ shows these precursor peaks in the repulsive Luttinger 
liquid for different system sizes. As the system size is 
increased the precursor peaks move towards lo — 0. For 
the calculation of the Drude weight these peaks should 
not contribute, and we use low-energy cut-offs to ignore 
them. 

Fig. |2l|and Fig. E2show the conductivity in the super- 
fluid phase in the repulsive and attractive Luttinger liq- 
uid regions. In the Luttinger liquid the conductivity was 
predicted to increase with a power-law for small fre qujen. - 
cies, and decay exponentially for big frequenciesEiltJ. 
The conductivity in the attractive region shown in Fig. 
p2| is in good qualitative agreement with this. In the re- 
pulsive case (Fig. |l]) there are too few peaks to clearly 
identify this behavior. Bigger systems would have to be 
studied to determine if the overall shape is qualitatively 
different from the attractive region. 



^ 2 - 




FIG. 21. The conductivity al^^{uj) in the repulsive re- 
gion of the superfluid phase at density p — 3/4, t = 0.05, 
U = 1, V = 0.4. Data shown is for different system sizes with 
(L = 64,m = 256,?? = 0.05), (L = 128,m = 256,r? = 0.05) and 
(i = 256,m = 512,r? = 0.2). Broadening rjg = 0.05 for the 
plot, and the low-frequency cut-offs given in Tab. IV 



The Drude weight in the repulsive and attractive 
regime of the superfluid phase is shown in Tab. H^ and 
Tab. |V[ For some system sizes data with different num- 
bers of states m and broadening 77 is shown. The numer- 
ical accuracy depends on these parameters, with bigger 
m and smaller 77 for higher accuracy. The data in Tab. 
IV and Tab. ^ shows that the impact of m and r] on the 
Drude weight is small. In both the attractive and repul- 
sive case we find big non-zero values that are close to the 
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FIG. 20. The precursor peaks in the conductivity al^^{uj) 
in the repulsive region of the superfluid phase at density 
p = 3/4, t = 0.05, U = 1, V = 0.4, system size L = 128 
and m — 128 states. Broadening rj = 0.02 for correction 
vectors, and rjg — 0.005 for the plot. 




FIG. 22. The conductivity a[''^{ij) in the attractive re- 
gion of the superfluid phase at density p — 3/4, t = 0.5, 
U = 1, V = 0.4. System size L = 32 and L = 64 with 
m — 128 states, and L = 128 with m = 256 states. Broaden- 
ing 77 — 0.2 for the correction vectors, rjg = 0.2 for the plot, 
and the low-frequency cut-offs given in Tab. M. 
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TABLE V. The Drude weight D for different system sizes 
in the attractive region of the superfluid phase at density 
p = 3/4, t — 0.5, U = 1,V — 0.4. Same notation as in 
Tab.Fn 



L D/t {T)/{tL) 2/tjol^'{Lo)du r, 



32 


1.438 


1.444 


0.006 


0.2 


128 


10"* 0.6 


64 


1.421 


1.427 


0.006 


0.2 


128 


10-* 0.4 


128 


1.412 


1.417 


0.005 


0.2 


128 


10'* 0.3 


128 


1.411 


1.417 


0.006 


0.2 


256 


10-^ 0.3 



kinetic energy per site in the systems. The differences 
in the Drude weight in different system sizes, with the 
exception of the smaUest systems, are rather due to nu- 
merical errors that grow with the system size, than due 
to finite-size effects. 



conductivity has the expected shape. In the repulsive 
region of the Luttinger liquid we found a different shape, 
but could not determine if this is due to the finite system 
sizes. The Drude weight or superfluid stiffness was found 
to be big in both the attractive and the repulsive region. 
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IX. CONCLUSIONS 



APPENDIX A: TRUNCATION OF THE MODEL 



In summary, we have studied the phase diagram of 
the one-dimensional Bose-Hubbard model with on-site 
only interactions and with additional nearest-neighbor 
interaction. The density matrix renormalization group 
(DMRG) was used to calculate chemical potentials for 
given densities and model parameters, and by doing this 
for sets of parameters the phase boundaries of the Mott- 
insulators and the charge density wave phase were deter- 
mined. 

The low-energy behavior of the superfluid phase of one- 
dimensional bosonic systems is that of a Luttinger liquid. 
We determined the Luttinger liquid parameter K from 
the decay of the hopping correlation functions. Since the 
value of K is known for insulator-supcrfluid transitions, 
we could use it to locate the Kosterlitz-Thouless transi- 
tions at the tips of the Mott-insulators and the charge 
density wave phase in the /i-t phase diagram. 

In the charge density wave phase we found that close 
to the phase transition the structure factor depends on a 
power-law of the superfluid correlation length. From this 
we conclude that there is a direct phase transition from 
the charge density wave phase to the supersolid, and no 
intermediate phase like a supersolid or normal phase. 

The charge density wave phase is surrounded by a re- 
gion of the superfluid phase where K > \, which cor- 
responds to a Luttinger liquid with repulsive effective 
interactions. Kane and Fisher have shown that this re- 
gion will be turned into an insulator by a single impurity. 
We determined the boundary of the repulsive region by 
finding the line where _R' = 1 in the phase diagram. We 
found that this boundary does not go to i = as the 
Mott-insulator is approached, but ends in the side of the 
Mott-insulating region, where the Luttinger liquid pa- 
rameter also in K — 1. 

We calculated the ac-conductivity in the Mott- 
insulator and the superfluid phase. In the Mott-insulator 
and the attractive region of the superfluid phase the ac- 



The number of possible states per site in the Bose- 
Hubbard model is infinite since there can be any number 
of particles on a site. For practical DMRG calculations 
we truncate the model by only allowing a maximuHji num- 
ber of particles Umax on each site. Pai et al.EZl chose 
nmax = '^lin El DMRG study, while Kashurnikov and 
SvistunovtJ used rimax = 3 in a Quantum Monte Carlo 
study. To verify the effect of this truncation on the corre- 
lation function V{r) = (blbo), we calculate systems with 
different Umax- Fig- |2^ shows T{r) in the Mott-insulator. 
Due to the small particle-hole excitations, the correlation 
functions are almost identical for Umax > 3. In the su- 
perfluid phase there are more particle-hole excitations 
which are affected by the truncation. Fig. e3 shows 
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FIG. 23. The r(r) = (b,.^o) correlation function for vari- 
ous truncations of the maximum number of particles per site 
rimax in the Mott-insulator at density p = 1 in a L = 128 
system with i = 0.1 and V = 0.4. For Umax > 3 the different 
correlation functions become indistinguishable. 
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FIG. 24. The r(r) — (blbo) correlation function for var- 
ious Umax in the superfluid phase in a L = 128 system with 
p = 1, t = 0.5 and V = 0.4. 



FIG. 26. The r(r) = (blbo) correlation function in the su- 
perfluid phase at p = 1 for various truncation errors. System 
size L = 128, t = 0.5 and V = 0.4. 



that the correlation functions are independent of nmax 
for Umax > 4. By choosing Umax = 5 for all calculations 
the effect of the truncation should be small enough not 
to affect the results presented in this work. 



APPENDIX B: TRUNCATION OF THE DMRG 
BASIS 

In every DMRG step the basis is truncated, and only 
the eigenstates of the density matrix with the biggest 
eigenvalues are kept. The density matrix weight of the 
discarded states A is a measure of the error caused by 
these truncations. To verify to which extent the trunca- 
tion errors affect the results, we calculate the correlation 



function T{r) — (blbo) with different numbers of states 
kept in the DMRG basis. Fig. H shows r(r) with dif- 
ferent truncation errors in the Mott-insulator. Even for 
very small numbers of states the discarded weight is very 
small, and the dependence on the weight of the discarded 
states is weak. Note that the discrepancies are mostly 
apparent due to the logarithmic scale. 

The correlation function in the superfluid phase is 
shown in Fig. g6[ We find that the discarded weight 
with the same number of states is bigger than it is in the 
insulator. At small distances r the correlation functions 
are very similar for all numbers of states, with increasing 
differences as r is increased. If the discarded weights A 
are smaller than 8 x 10~^, the correlation functions coin- 
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FIG. 25. The T{r) — (blbo) correlation function for var- 
ious truncation errors A in the p — 1 Mott-insulator phase. 
System size L = 128, t = 0.1 and V = 0.4. 
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FIG. 27. fi versus the discarded weight A in the 

Mott-insulator at p = 1. The scale for A is logarithmic on 
the left plot and linear on the right plot. The system size is 
L = 128, t = 0.1 and V = 0.4. 
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FIG. 28. fi'^ and fi~ versus discarded weight A in the 
superfluid phase. The scale for A is logarithmic on the left 
plot and linear on the right plot. The density is p = 1, system 
size L = 128 and t = 0.5 and V = 0.4. 



cide even at the boundaries of the system. By requiring 
the discarded weight to be smaller than A < 10~^ for the 
calculation of the correlation functions, accuracy should 
be high enough in all cases. 

The chemical potentials are calculated from the ener- 
gies it takes to add a particle or hole. Figures |^ and |28| 
show chemical potentials calculated with different num- 
bers of states kept versus the discarded weight A. The er- 
ror bars correspond to the changes in the energies during 
a DMRG sweep. Differences in the chemical potentials 
are small for A < 10~^. We require the discarded weight 
to be smaller than A < 5 x 10~^ for the calculations of 
the chemical potentials, and to improve the results fur- 
ther, we extrapolate linearly from the two values with 
the lowest A. 



APPENDIX C: STRONG COUPLING 

EXPANSION AT THE 

COMMENSURATE-INCOMMENSURATE 

TRANSITION 

The fundamental difference between 

the commensurate-incommensurate phase transition at 
p = 1/2 in one and two dimensions can be illustrated 
with the help of a strong coupling expansion. In the 
strong coupling limit the kinetic energy is zero. The zero 
order states are the groundstates of the Hamiltonian only 
including the particle-particle repulsion: 



H = U^ni{ni - l)/2 + V^niUi+i 



(CI) 



FIG. 29. The CDW in 2-d at t = 0. 



H' ^ -tJ2ib]b^ + h.c.) . 



(C2) 



Strong coupling expansions of this type have been suc- 
cessfully used tp-staidy the phase diagram with on-site 
only interactionEjoa. To determine the phase bound- 
aries of the Mott-insulator at p = 1, first the groundstate 
of Eq. ( |Cl| ) has to be found. In this state there is simply 
one boson sitting on every site. Higher terms of the per- 
turbation series introduce local particle-hole excitations. 
The chemical potentials on the boundaries can be deter- 
mined from the energy it costs to add a particle (Eq. (||)) 
or a hole (Eq. (0)). But in these cases, the zeroth order 
ground state is degenerate, since the additional particle 
or hole can sit on any site. This degeneracy is lifted in 
first order perturbation theory. In first order the prob- 
lem is reduced to the additional single particle moving 
on a uniform background of completely localized parti- 
cles. Since the extra particle gains energy by hopping 
from site to site, it becomes completely delocalized. Al- 
though this behavior can be modified in higher order of 
the perturbation series, and there are limitations due to 
the radius of convergence, it is interesting to note the 
difference between the perturbation series in the insula- 
tor and with an additional particle or hole. At integer 
density the series starts out with a completely localized 
state, while it starts with a completely delocalized state 
if there is an additional particle or hole. This is in good 
agreement with the fact that there is a Mott-insulator for 
small t at integer density, and a direct phase transition to 
a superfluid (delocalized) phase if the density is changed. 

A similar strong coupling expansion can also be used 
at the charge density wave phase at p = 1/2. The charge 
density wave phase at t = is a state with alternating 
particle numbers, in one as well as two dimensions (Fig.|29| 
and Fig.pffl). Higher order terms in the perturbation se- 
ries introduce local particle hopping without destroying 
the charge density wave order. 



The series expansion is made in terms of the kinetic en- 
ergy term: 



FIG. 30. The CDW in 1-d at i = 0. 
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FIG. 31. CDW with additional particle in 2-d at f = 0. 



At the comniensurate-incomniensurate transition, an 
additional particle (or hole) enters the system. For 
V < U/2 the energy is smallest if the additional particle 
goes to one of the empty sites. Fig. ^ shows how the ad- 
ditional particle fits into the two-dimensional charge den- 
sity wave. In higher orders of perturbation theory the ad- 
ditional particle, as well as particle-hole excitations, hop 
on the charge-density background without destroying it. 
From this we cannot infer if the true (non-perturbation 
theory) groundstate is superfluid or not, and if the charge 
density wave order is destroyed by the particle hopping. 
Nevertheless, it is interestipg-tp note that for this case su- 
persolids have been foundEJO in two dimensions. Close 
to the charge density wave phase at p = 1/2, the charge 
density order survives at small doping. 

In contrast to this, the one-dimensional case looks 
quite different. The additional particle also goes to an 
unoccupied site. If the charge density wave remains un- 
changed, the additional energy is AE = 2V. With the 
structure factor S^r, the charge density wave is given by 
an order parameter (ni) = p + S*^ exp (jtt^ -I- i(/<o)- An 
additional particle or hole can also be added by shifting 
the phase do by n over a region with an odd number of 
sites. Fig.g^ shows an example of such a state. In the 
center there is a domain with a tt phase shift, and the 
number of particles compared to the charge density wave 
(Fig.ES) is increased by one. And the additional energy is 
also AE = 2V. To lift the degeneracy between all these 
states in first order perturbation theory, it can be seen 
that a particle hopping next to the domain boundary is 
equivalent to the domain wall moving. Since energy is 
gained by this, the domain walls are completely delocal- 
ized in first order perturbation theory. Unlike the two 
dimensional case, where the charge density wave order 
survives in all orders, in one dimension it is destroyed in 
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FIG. 32. Additional particle in 1-d at i = 0. The thick 
lines are the domain walls, i^o is the phase of the charge 
density wave. 



first order. While a perturbation series does not neces- 
sarily converge, this striking feature illustrates the funda- 
mental difference between the one and two dimensional 
case. 



APPENDIX D: CURRENT OPERATOR WITH 
OPEN BOUNDARIES 

DMRG calculations work best with open boundary 
conditions. To calculate the conductivity with DMRG, 
the current operator has to be implemented. The current 
operator as it is given in Eq. ( p^ ) can be used directly 
with periodic boundary conditions, but to apply the cur- 
rent operator with open boundary conditions we modify 
it with a filter function: 



Jq=0 



E 

n— — oo 



P{x,JM){bl.^b„-h.c.) 



(Dl) 



The filter function P{xn/M) used here is a Parzen filter, 
Xn is the distance of site n from the middle of the system, 
and M — L/2 is half the system size. The Parzen filter 
looks very similar to a Gauss function, but goes smoothly 
to zero at the system boundaries. It is given as: 



P{x) 



l-6|a;|2 + 6|a;|3 if < \x\ < 1/2 
2(1- \x\f if 1/2 < |x| < 1 



(D2) 



A prefactor a is chosen to provide results with the same 
amplitude as those found in systems with periodic bound- 
ary conditions, with a chosen so that ^P(x„/M)^ = 1. 
To verify the effect of open boundaries, we do two sep- 
arate calculations of the conductivity, one with periodic 
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FIG. 33. The conductivity (t['^®(cij) in the insulator with 
periodic and open boundary conditions. The density is p = 1, 
t — 0.05, [/ = 1 and V = 0.4, system size L — 32, and 
771 — 128 states. Broadening rj — 0.2 for the correction vec- 
tors, and Tjg = 0.01 for the plot. The inset shows the data on 
an expanded scale. 
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FIG. 34. The conductivity al''^{u;) in the superfluid 
with periodic and open boundary conditions. The density 
is p = 3/4, t = 0.05, U = 1 and V = 0.4, system size L = 32, 
and m = 128 states. Broadening rj — 0.2 for the correction 
vectors, and rjg = 0.02 for the plot. The inset shows the data 
on an expanded scale. 



and one with open boundaries, for otherwise identical 
system parameters. Fig. B3^ and Fig. M show the con- 
ductivity in the Mott-insulator and the superfluid phase 
with open and with periodic boundary conditions. Even 
in the small systems with L = 32 the curves are quite 
similar. In the superfluid phase a precursor peak at 
small frequencies can be seen in the system with open 
boundary conditions. Fig. pfl shows how the precur- 
sor peaks move to smaller frequencies as the system size 
is increased. They are an artifact of the open bound- 
ary conditions, and we use frequency cut-offs to exclude 



them from the calculation of the Drude weight. Tab. VI 
shows the Drude weights in the insulator and the super- 
fluid. The values for open and periodic boundary con- 
ditions compare quite well, and we estimate an error of 
AD/t « 0.02. 

TABLE VI. The Drude weight D with open and periodic 
boundary conditions. Systems size L — 32, t = 0.05, U = 1 
and V — 0.4, m — 128 states, broadening rj — 0.2. Same 
notation as in Tab.[lll[ 

boundary p D/t {T)/(tL) 2/t J cr[''^{uj)duj A ujc 



periodic 1 0.0195 0.6392 0.6197 

open 1 0.0017 0.6392 0.6375 

periodic 3/4 0.8050 0.8957 0.0907 

open 3/4 0.7881 0.8836 0.0955 
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